library(tidyverse)
library(sf)
library(here)

# Function to convert degrees to radians
deg2rad <- function(deg) {
  return(deg * pi / 180)
}

# Unit normal function in R
unit_normal <- function(elevation_deg, azimuth_deg) {
  elev_rad <- deg2rad(elevation_deg)
  azim_rad <- deg2rad(azimuth_deg)
  
  return(c(sin(azim_rad) * sin(elev_rad),
           cos(azim_rad) * sin(elev_rad),
           cos(elev_rad)))
}

# Alpha_s function in R
alpha_s <- function(alpha_i, direct_fraction, s_hat, n_hat, n_elevation_deg) {
  num <- (1 + direct_fraction * (s_hat[3] - 1))  # s_hat[3] = z component in R
  denom <- (direct_fraction * max(sum(n_hat * s_hat), 0) +
              (1 - direct_fraction) * (180 - n_elevation_deg) / 180)
  
  if (denom == 0.0) {
    return(alpha_i)
  }
  return(alpha_i * num / denom)
}

# Compute segment surface albedos function in R
compute_segment_surface_albedos <- function(solarZenithDeg, solarAzimuthDeg, albedoMedian, normal, pitchDegrees) {
  s_hat <- unit_normal(solarZenithDeg, solarAzimuthDeg)
  
  alpha_s = alpha_s(alpha_i = albedoMedian,
                    direct_fraction = 0.85,
                    s_hat = s_hat,
                    n_hat = normal,
                    n_elevation_deg = pitchDegrees)
  
  return(alpha_s)
}

# surface_albedo_test <- roof_discs %>% 
#   drop_na(azimuthDegrees) %>% 
#   sample_n(20) %>% 
#   select(NEON_solar_elevation_mode, NEON_solar_azimuth_mode, NEON_albedo_median, pitchDegrees, azimuthDegrees) %>% 
#   rowwise() %>% 
#   mutate(normal = list(unit_normal(pitchDegrees, azimuthDegrees)),
#          NEON_alpha_s = compute_segment_surface_albedos(solarZenithDeg = 90 - NEON_solar_elevation_mode, 
#                                                         solarAzimuthDeg = NEON_solar_azimuth_mode, 
#                                                         albedoMedian = NEON_albedo_median, 
#                                                         normal = normal, 
#                                                         pitchDegrees = pitchDegrees))
# 
# write_csv(select(surface_albedo_test, !normal), here("data", "surface-albedo", "surface-albedo-test.csv"))
# python_test <- read_csv(here("data", "surface-albedo", "python-surface-albedo-test.csv")) 
# 
# test <- surface_albedo_test %>% 
#   select(NEON_alpha_s) %>% 
#   bind_cols(select(python_test, NEON_alpha_s))

# test_discs <- st_read(here("data", "paper-analysis", "boulder", "roof-test-discs.geojson")) %>% 
#   select(ID = id) %>% 
#   filter(ID == "00000000000000005041") %>% 
#   left_join(roof_discs_alb)
# 
# x <- test_discs %>% filter(PNeo_date == "2024-05-16_17-58-54" & albedo_date == "2024-05-16")
# 
# n_hat <- unit_normal(x$pitchDegrees, x$azimuthDegrees)
# print(n_hat)
# 
# s_hat <- unit_normal(90 - x$solarElevation, x$solarAzimuth)
# print(s_hat)
# 
# test <- compute_segment_surface_albedos(90 - x$solarElevation,
#                                         x$solarAzimuth,
#                                         x$Med_albedo,
#                                         n_hat,
#                                         x$pitchDegrees)
# print(test)
